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APPLICATION OF BOUNDARY INTEGRAL METHOD TO ELASTOPLASTIC 
ANALYSIS OFV-NOTCHED BEAMS 


by Walter RzasnickI, Alexander Mendelson, Lynn U. Albers, 
and Demetrius D. Raftopoulos* 

Lewis Research Center 

SUMMARY 

The boundary integral equation method was applied in the solution of the plane elas- 
toplastic problem. The use of this method was illustrated by obtaining stress and strain 
distributions for a niunber of specimens with a single edge notch and subjected to pure 
bending. The boundary integral equation method reduced the inhomogeneous bihar- 
monic equation to two coupled Fredholm-type integral equations. These integral equa- 
tions were replaced by a S 3 ^tem of simultaneous algebraic equations and solved numer- 
ically in conjunction with a method of successive elastic solutions. 


INTRODUCTION 

Knowledge of the stress distribution in the neighborhood of a singularity, such as the 
tip of a V-notch in a bar loaded in tension or bending, is of fundamental importance in 
evaluating the resistance to fracture of structural materials. Elastic solutions to vari- 
ous geometries have been obtained by a number of different methods. Among the more 
effective ones, are the complex variable method (ref. 1), collocation method (ref. 2), 
and finite element method (ref. 3). However, the first two of these methods are not gen- 
eral enough nor readily adaptable to three-dimensional or elastoplastic problems. And 
the finite element method requires solutions of large sets of equations and fails to give 
sufficiently fine resolution in the vicinity of notch tips. 

The recently developed boundary integral methods (ref. 4) offer an attractive alter- 
native to other methods of analysis. These methods have a number of advantages as 
listed in references 4 and 5, the most important of which is that nodal points are needed 
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only on the boundary instead of throughout the region as required by finite element 
methods. To date these methods have been used primarily for obtaining elastic solu- 
tions to various problems. Their extension to elastoplastic problems has been proposed 
in references 4 and 6. However, no elastoplastic solution has heretofore been obtained. 

Reference 5 describes in detail the application of a boimdary integral method to the 
solution of the elastic problem of a V-notched beam in pure bending. The present re- 
port extends this method to the more complicated elastoplastic problems. Solutions of 
such problems by finite elements has not been too successful in obtaining fine enough 
resolution and sufficiently accurate results in the vicinity of the notch tip (ref. 7). The 
present method overcomes these difficulties, since the stress and strain can be com- 
puted at any arbitrary point in the body from a knowledge of boundary values only'. 
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SYMBOLS 

area of cell in region R 
coefficients in stress equations (14) 

notch depth 

coefficients inboundary equations (10) 

dimensionless notch depth, a/w 

boundary contour 

Young’s modulus of elasticity 

function of plastic strain increments 

path independent contour integral 

stress intensity factor for mode I 

generalized stress intensity factor for mode I 

half length of beam 

strain hardening parameter, ratio of slope of linear strain-hardening 
curve to elastic modulus 

order of stress singularity 

unit vectors normal and tangent to contour C 
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point on contour C or in region R 
point on contour C 
dimensionless load, a„„ /a 

XliaX (J 

planar region bounded by closed contour C 

distance between two points having coordinates (x, y)^ and U>T?)j 

polar coordinates 

length measured along contour C 

convergence parameter 

displacement vector 

displacement in y direction 

width of beam 

rectangular Cartesian coordinates 

dimensionless rectangular Cartesian coordinates, x/w, y/w 
notch angle 

^2 g2 

Laplace's operator, + 
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biharmonic operator, + 2 + 
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Kronecker delta 
strain tensor 

equivalent modified total strain 

equivalent plastic strain increment 

components of strain tensor in Cartesian coordinates 

modified total strain tensor 

components of modified total strain tensor in Cartesian coordinates 
Poisson's ratio 

rectangular Cartesian coordinates 
2 

function of r, r In r 

equivalent stress 

maximum nominal bending stress 
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components of stress tensor in Cartesian coordinates 
tensile yield stress 

2 

function of Airy stress function, V <p 
Airy stress function 


Subscript: 

i, j,k integers 

Superscripts: 
e elastic 

p plastic 

~ dimensionless quantity 

' derivative in outward normal direction, 9/9n 


METHOD OF SOLUTION 

The problem of determining the state of stress and strain in a plane elastoplastic 
problem can be reduced to solving the following inhomogeneous biharmonic equation 
for the Airy stress function, as shown in reference 8: 


V^(p = g(x, y) 
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for the plane strain case and 


g(x, y) = -E 
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for the plane stress case, where, eL eP, andeP represent the accumulation of plastic 

X y xy 

strain increments from the beginning of the loading history up to, but not including the 
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current increment of the load, and and ^e^y are the increments of plastic 

strain due to the current increment of load. 

The stress function cp must satisfy appropriate boundary conditions. For the 
problem under consideration (fig. 1), <p(x, y) and its outward normal derivative dcp/dn 
must satisfy the following boundary conditions (ref. 9): 


<p(x, y) = 0; -^ = 0 along boundary OA and OA* 
3n 


^(x, y) = 0; -^=0 aloi^ boundary AB and A'B' 
3n 






(4) 


■^ = 0 along boundary BC and B’C' 
0n 


a w 


(p{x, y) = — ; .^ = 0 along boundary CD and C’D 

6 0n 


To solve equation (1) by means of the boundary integral method, use is made of 
Green's second theorem to reduce this equation to coupled integral equations, as shown 
in references 4 and 9. The result is 
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and r(x, y; tj) is the distance between any two points P(x, y) and q(^, rj) in the re- 
gion R bounded by the curve C, such that P C R + C and q C C (fig. 2). 

Equation (5) would, for a known function g(x, y), give us directly a solution to the 

p 

biharmonic equation (1) provided the functions ^(x, y), d(p{x,y)/dn, V^<^(x, y), and 

a[v2</>(x, y)]/9n were known on the boundary C. 

However, only the stress function <p and its outward normal derivative dcp/ 9n are 

2 2 

specified (eq. (4)). The values of V (p = ^ and 8(V ^)/9n = 04>/0n on the boundary 
must be compatible with the given values of (p and dcp/dn. To assure this compatibil- 
ity, we have to solve the system of coupled integral equations (6) and (8), which contain 
the unknown functions $ and 04>/0n. 

Once the values of $ and 04>/0n on the boundary C of region R are known we 
can proceed with the calculation of the stress field in the region R utilizing equation (5) 
and the equations which define (p, namely. 
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The calculation of the function g(x, y), which is obtained iteratively, will be dis- 
cussed subsequently. 


NUMERICAL PROCEDURES 
Solution of the Integral Equations 

Since it is generally impossible to solve the system of coupled integral equations 
analytically, a numerical method is utilized in which the integral equations (6) and (8) 
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are replaced by a system of simultaneous algebraic equations. 

For simplicity of notation the normal derivatives are denoted by prime superscripts. 
The boundary is divided into n intervals, not necessarily equal, numbered consecu- 
tively in the direction of increasing s. The center of each interval is designated as a 
node. The values of and are assumed constant on each interval and equal to the 
values calculated at the node. 

In similar manner the interior of region R is covered by a grid, containing m 
cells. The cells do not have to have equal areas. Their nodal points are located at the 
centroids. The value of g(^, ?]) is assumed constant over each cell and equal to the value 
calculated at the centroid. The arrangement of boundary and interior subdivisions is 
shown in figures 3 and 4. 

Using these assumptions, equations (6) and (8) can be replaced by a system of 2n 
simultaneous algebraic equations with 2n unknowns, that is, and 
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where integration is taken over the j interval, and r.. is the distance from i^" node 

th ^ 

to any point in the j interval. The normal derivatives in equations (11) are taken on 
the interval. 

For curved boundaries the coefficients given by equations (11) can be evaluated, if 
necessary, by Simpson's rule for i j. For i = j, because of the singular nature of 
the integrand, the integrals for the coefficients must be evaluated by a limiting process. 
For boundary intervals, such as for the problem treated herein, which can be repre- 
sented by straight lines a closed form solution can be obtained for these coefficients. 
Boundary equations (10) expressed in matrix form become 
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(12) 


Thus, the problem is reduced to the solution of the following matrix system; 


[B]{X} = {R} (13) 

where [B] is 2n x 2n matrix and {x} and {R} are 2n x 1 column matrices. 

Matrix [b] is dependent only on geometry, that is, number of nodes and their dis- 
tribution on the boundary. Since the matrix (r) contains the nonlinear function g(^, r\), 
which depends on the stress field and therefore on matrix (x}, an iterative process will 
be used to obtain the solution. 

To calculate stresses, at any nodal point in the region R, from the stress function 
cp, we need not perform any numerical differentiation. Equation (5) can be differentiated 
under the integral sign and once $ and are known on the boundary the stresses can 
be obtained by the same type of numerical integration as in equations (10). Applying 
equations (9) to equation (5) yields for the case of a rectangular grid the following stress 
equations: 
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8;rcj^(x, y). = 



where now i = 1, 2, 3, . . . , m and refers to the centroid of the cell, 6^^ and 6y 
represent, respectively, x-directional and y-directional dimension of the cell. The 
coefficients A^j, B^, C-j, D^, E.., F^, G^, H^, I.., and K-j are obtained by appro- 
priate differentiation under the integral sign of the coefficients given by equations (11) 
and are listed in appendix A. 

The stress function is not constant on the loaded boundaries BC and B'C’. The 
assumption that it is piece-wise constant may lead to appreciable errors in the numerical 
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results. To eliminate this source of error, the summations given in equations (10) 
and (14) for intervals lying on the loaded boundaries and involving the stress function are 
replaced by direct integration. 


Boundary Interval and Interior Grid Size 

The number of nodal points prescribed for the boundary is theoretically unlimited. 
However, computer storage capacity for the computer used and difficulties associated 
with inversion of large matrices limited the order of the coefficient matrix [b] of equa- 
tion (13) used herein to 140. 

Because of geometric and loadii^ symmetry about the x-axis, it is possible to re- 
duce the number of unknowns. For 2n total number of nodal points the number of equa- 
tions and unknowns, and $!, is reduced from 4n to 2n. Additional reduction in the 
number of unknowns is accomplished by taking into consideration St. Venant's effect at 
the loaded boundaries (ref. 2). 

Since the vicinity of the crack tip is of greatest interest, a fine nodal spacing along 
the notch was chosen. To reduce the error introduced by the change in the interval size 
(ref. 10) around boundary points A and A' and at the same time to obtain fine resolu- 
tion at the tip of the notch, the boundary along the notch was divided into a number of 
intervals progressively decreasing in length. The rate of change in the interval length 
and the resulting ler^h of the smallest interval was found to have a great influence on 
the stress field in the vicinity of the tip of the notch. The rate of change in the interval's 
length alor^ these boundaries was optimized by the method presented in reference! 5. 

For the cases considered optimum ratios of the lengths of two consecutive boundary in- 
tervals were found to be in the range of 1. 08 to 1. 10. The resulting smallest dimension- 
less boundary interval length varied from 0. 0001 to 0. 0002. The nodal arrangement 
shown in figure 5 was used for all cases considered, resulting in a set of 140 equations 
containing 140 unknowns. Note that the corner points are always designated as interval 
points, never as nodal points, thus eliminating discontinuous functions from munerical 
analysis. 

The choice of the size of the grid, which has to cover the region where plastic flow 
is expected to occur, is of utmost importance. A too coarse grid will not detect changes 
in the values of plastic -strain for small loading increments. A too fine mesh size may 
result in distorted values of second-order derivatives of plastic strains, which appear 
in the function g(x, y). The loading increment and the grid size are related to each 
other. A bad choice of either of them could result in the divergence of the iterative 
process. To allow the maximum of grid points to be within the expected plastic zone, a 
variable grid spacing was chosen. The grid used for plane strain conditions was finer, 
in general, than the one used for plane stress case. 
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The interior region, where plastic flow is expected, was divided into r x s rec- 
tangular cells. Due to symmetry about the x-axis, the number of unknown functions g, 
appearing in the boundary equations (10) and stress equations (14), was reduced from 
rxs to m = rx(s + l)/2, where now the coefficients of these functions represent the 
sum of the effect of left-hand and right-hand sides of the plastic field. Because of com- 
putation time limitations, the grid was arranged in 27 x 23 cell pattern, resulting in 
the number of unknowns g to be equal to 324. By increasing the number of unknowns to 
400, the computation time for one iteration almost doubled. The smallest cells, located 

in the vicinity of the tip of the notch, have dimensions = 0. 004; 6,,/w = 0. 008 for 

X y 

plane strain cases, and 6.^,/w = 0. 004, 6,,/w = 0. 016 for plane stress cases. 

X y 

Cells with centroids outside the boundary of the beam above the tip of the notch were 
discarded and corresponding values of g were set equal to zero. A typical interior grid 
is shown in figure 5. 

Central difference equations, given in appendix B, were used to evaluate g(x, y) 
only for interior plastic cells, that is, where there was plastic flow at all eight adjacent 
cells. For noninterior cells the function was taken as an average of values of g at 
neighboring plastic cells. Other methods of dealing with g at centroids of noninterior 
plastic cells, such as backward differences of extrapolation, led to oscillation or diver- 
gences of the iterative process. 


Method of Successive Elastic Solutions 

The solution to this boundary value problem is obtained by an iterative process, 
known as the method of successive elastic solutions described in detail in reference 8. 
This method, applied to the present problem, proceeds as follows. The loading path is 
divided into a number of sufficiently small increments Act .. The g function in the 

HlaX^ 1 

inhomogeneous equation (1) is thought of as a sum of Ag^ functions, each corresponding 
to its load increment Act Each Ag. is defined in terms of derivatives of the cur- 

rent plastic strain increments, Ae^, Ae^, and Ae^y, which in turn depend, through 
equivalent plastic strain increment Ae^, on the equivalent stress ct^ associated 
with the previous load. The current plastic strain increments change iteratively as 
changing Ag^ affects the stress field. 

The iterative procedure for determining plastic strain increments for each load in- 
crement is as follows: 

(1) Select a value of load cr .. 

IXloX^ 1 

(2) Guess initial values of plastic strain increments. For the first load increment, 
assume all values to be equal to zero. Otherwise, use the converged values from the 
previous load increment. 

(3) Calculate g^ = gi_i + Agp 
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(4) Calculate 4*^ and 4>! from boundary equations (10). 

(5) Calculate the stresses from stress equations (14). 

(6) Calculate modified total strains e|j and the equivalent modified total strain 
from the following equations: 


» e * D 
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where are the elastic strains computed from the stresses obtained in step (5). 

(7) Calculate the equivalent plastic strain increment from 
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for plane stress case, where ^ is the value of the equivalent stress at the end of 
the previous increment of loading. For the first load increment and also for the case 
where there was no plastic flow under previous loading, ._j is equal to the yield 
stress Oq. For the case of linear strain hardening equation (17) becomes 
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(8) Calculate a new set of plastic strain increments from 
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(9) Repeat steps (3) to (8) imtil the plastic strain increments converge. 

(10) Svun the plastic strain increments and return to step (1). 

It should be noted at this point that the aforementioned procedure can be applied only 
where there is no unloading. Once the successive approximation procedure has con- 
verged, the stresses and strains are known everywhere in the beam. 

The iterative process is illustrated by the flow diagram of figure 6. 


Convergence 


The convergence criterion is defined as an arbitrary maximum difference between 
successive values of one or more iterants. For the plastic field containing many points 
at which convergence is required, it is reasonable to set the convergence criterion 
based on an average value of the change in plastic strain increments. For the problems 
considered in this report, the convergence criterion was based on the convergence of 
plastic strain increments Ae^ and was defined as 
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where k refers to current iteration, n refers to number of plastic grid points for cur- 
rent iteration, and T is an arbitrary convergence parameter. 
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The computations were performed on a digital computer using a FORTRAN IV pro- 
gram with single-precision arithmetic. The matrix system given by equation (13) was 
solved using the modified Gauss elimination method, which utilizes pivoting and forward 
and backward substitutions. 

The convergence parameter T was set to 0. 003 or 0. 005 depending on the load in- 
crement. Decreasii^ the convergence parameter T to 0. 001 resulted in change in 
plastic strain increment values of approximately one in the third or fourth significant 
figure. For higher order accuracy, which does not appear to be necessary from a prac- 
tical point of view, the computation time would be prohibitively long. The fact that the 
method of successive elastic solutions converges to the right answer has been shown by 
many examples in references 8 and 11. 


RESULTS AND DISCUSSION 

A number of beam problems were solved for both plane strain and plane stress 
cases. These included notch depth to beam depth ratios of 0. 3 and 0. 5, notch armies of 
3° and 10°, strain hardening parameter values of 0. 05 and 0. 10. In addition, calcula- 
tions were performed using the actual stress-strain curve of a 5083-0 aluminum alloy. 
For all cases Poisson’s ratio was set at 0. 33. 

The load increment size used was necessarily a compromise between the accuracy 
desired and computational time required for convergence. For strain hardenir^ param- 
eter m = 0. 05 the load increment size Aq was taken equal to 0. 05; while for m = 0. 10, 
Aq = 0. 10. For the case of a 5083-0 aluminum alloy, where the actual stress-strain 
curve (fig. 25) was used, the load was incremented by Aq = 0. 025. 

For the beam with dimensionless notch depth a = 0. 5 the minimum load required 
to produce plastic flow at the most highly stressed grid points was found to be q = 0. 30, 
and for a = 0. 3 the initial load was found to be q = 0. 50. The maximum load consid- 
ered was q = 0. 7 for the a = 0. 5 cases, and q = 0.9 for the a = 0. 3 cases. In the 
process of solving the aforementioned problems, the case with strain hardening param- 
eter m = 0. 05 required approximately 50 iterations for each increment of load (i. e. , 

Aq = 0. 05) for the relatively fine convergence parameter used. For cases where the 
strain hardening parameter m = 0. 10 the average number of iterations needed for each 
increment of load (i. e. , A q = 0. 10) was reduced to 40, while use of the actual stress- 
strain curve resulted in convergence in approximately 10 iterations for the plane strain 
case and in 20 iterations for the plane stress case. 

Typical results of the computations are presented in figures 7 to 27 and tables I 
and n. Complete detailed results are given in reference 9. 

The growth of the plastic zone with load is shown in figures 7 to 14. It is seen that 
the shapes of the elastoplastic boundaries remain similar to each other as the load in- 
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creases. As expected, plastic flow starts around the tip of the notch and as the load in- 
creases appears also at the boundary opposite the notch. Comparison of figures 11 
and 12 with figures 13 and 14 shows that for the same loads the size of the plastic zones 
for plane strain are considerably smaller than for plane stress. 

The equivalent stress contours in the vicinity of the notch for maximum applied 
loads are plotted for two cases in figures 15 and 16. The curves are the loci of all 
points of constant equivalent stress. The curves corresponding to CTg/o^O “ ^ indicate 
the boundary of the plastic zone. In addition, an elastic yield locus representing the 
elastoplastic boundary based on the elastic solution is shown in each case. Since this 
is commonly assumed to be the boundary of the plastic zone, we can see that for plane 
strain cases this assumption introduces considerable error along the x-axis. Aloi^ 
this axis the lengths of plastic zones obtained by elastoplastic and elastic solutions vary 
by a factor of about two. 

Stresses and strains were calculated in all cases for interior grid points. Typical 
results of these calculations are given in figures 17 to 22. 

The order of the stress singularity n was determined for each case by the method 
described in reference 5. The results are given in tables I and n. In case of plane 
strain conditions, the stress singularity decreases slowly as loading increases. In the 
case of the plane stress condition we have a sudden drop in n from its elastic value as 
plastic flow appears. Subsequently n slowly increases approaching a limit as the load 
increases. 

The product of y-directional stress and total strain was calculated for various 
cases. The order of singularity of that product was determined by plotting In(ayey) 
against In r and by making a least squares fit of a straight line through the plotted 
points. It was found to be very close to unity for all cases considered. 

In the case of an elastoplastic problem the stress intensity factor Kj defined in 
reference 2, must be generalized to the form 


K?(a 


I' max 


) = limV^ r 
r-0 


c^y(r, e) 


9=0 


(23) 


For linear elastic behavior is identical with Kj. 

Variation of the dimensionless generalized stress intensity factor with load is shown 
in figure 23 for the case of a specimen with notch depth of a = 0. 5 and a = 10°, under 
plane strain condition and two values of strain hardening parameter m. The stress 
intensity factor shows no significant increase over the linear elastic value up to an ap- 
plied load of q = 0. 40. Above this load increases progressively for both m's, at 
the faster rate for lower strain hardening parameter. 
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The y-directional notch opening displacement was obtained for each case by numer- 
ical integration of the relation e.. = (l/2)(u. . + u. .) along straight line paths. For 
each case a number of paths were chosen through the plastic region near the notch, and 
the resulting displacements were averaged. In general, the notch opening displacement 
varies linearly with the load until the plastic zone is established at the boundary oppo- 
site the notch. Then it increases rapidly, reaching values several times that which 
would be calculated from the elastic solution. 

In order to verify in part the accuracy of the method used, a comparison of notch 
opening displacements was made with experimental results obtained by Bubsey and Jones 
(private communication from R. T. Bubsey and M. H. Jones of NASA Lewis Research 
Center). The specimen used in this experiment, made of aluminum 5083-0 with a length 
to width ratio of 4 and a crack length a = 0. 5, was subjected to three-point bending. The 
stress-strain curve for this specimen is shown in figure 25. The experimental results 
as shown in figure 24 are in good agreement with numerical results obtained herein. 

Finally, the J integral was evaluated for several cases by using relations given in 
reference 5. As in notch opening displacement calculations, straight line paths were 
chosen through the plastic zone near the tip of the notch. The integral was evaluated 
using values of stresses, strains, and displacements at cell centroids for a number of 
paths. The path independence of J was not conclusive, since the results varied up to 
15 percent from the averaged value. It is possible that the results obtained herein do 
not indicate that the path independent property is lost but rather that the field values of 
the displacements are not calculated with sufficient accuracy. 

The average values of the dimensionless J integral as a function of load are plotted 
in figure 26 for a case of a specimen with a 10° edge notch, a = 0. 5, m = 0. 05, and plane 
strain condition. At the start of plastic flow J^ increases rapidly with load. This is 
followed by almost linear variation with additional load. 

The relations between the J integral and stress intensity factor Kj developed for 
linear elasticity are obviously not applicable for the elastoplastic problem. By plotting 
the J/Kj ratios as a function of load q, the relation between the J integral and the 
dimensionless generalized stress intensity factor is obtained. Typical plots are 
shown in figure 27. In all cases, the ratio J/Kj^ remains almost equal to elastic value 
of 0. 89 for plane strain or 1. 0 for plane stress and increases sharply at the load corre- 
sponding to the appearance of the plastic zone at the boundary opposite the notch. Once 
the transition occurs the ratio increases approximately proportionally to the load incre- 
ment. 


16 



CONCLUSIONS 


The boundary integral equation method proved to be capable of giving very detailed 
results such as stress and strain distributions around the tip of the notch and, related 
to them, the shapes of plastic zones. This was accomplished using a relatively small 
number of unknowns. 

The obtained results also provide the information on the effect of strain hardening 
and on the differences that occur between plane stress and plane strain solutions. The 
singular nature of stresses and strains in the vicinity of the tip of the notch was con- 
firmed. The order of singularity for the strain energy density was found to be unity, 
which is consistent with results previously obtained by other investigators. 

The generalized stress intensity factor was introduced and calculated for several 
cases. The path independence property of the J integral was qualitatively confirmed 
and the relation between J and the generalized stress intensity factor was graphically 
extended to materials deforming according to the Prandtl-Reuss theory of plasticity. 

The presence of a singularity at the tip of the notch makes accurate answers very 
difficult to obtain. Nevertheless good agreement was obtained between the calculated 
results and experimentally measured notch opening displacement as shown in figure 24. 
Some improvement in the solution techniques and further investigation of the influence of 
the boundary nodal spacing and interior grid size on the resulting stress and strain 
fields, and therefore, on the notch opening displacements and J integrals, may be 
desirable. 

Lewis Research Center, 

National Aeronautics and Space Aclministration, 

Cleveland, Ohio, January 7, 1974, 

501-21. 
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APPENDIX A 


COEFFICIENTS OF THE STRESS EQUATIONS 


The coefficients appearing in stress equations (14) are given by the followii^ rela- 


A„=^(e„)=4 


■f 


9 (Xi - 1)2 - (y. - t?)2 


(x^ - ^)2 + (y^ - 77)2] 


32 / (yj - v) - (x^ - 

^ (yi - 4 




2(y4 - v) 

+ Vds 

(x. - ^)2 + (y. - r?)2^ 




lJi[(Xi - 4)2 + (y. - 77)2] 


2(yi - nr 1 

+ 1 + i^ds 

(x^ - 4)^ + (Yi - nf J 




2(x. - 4)' 


+ i Vds 

(x. - 4)2 + (y. - 77)2^ 


Ei.=^(d,.) = - 

dx^ 


- + (Yi - nf\ 


2(x. - 4)^ 

(Xi - 4)^ + (Yi - n)2 


+ l>ds 
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°ii = 


3x 3y 



ds 


H =^_(f ) = 8 
^ 0x 0y J 



(xj - |)(y^ - v) 


■ds 


( [(X^ - 1)^ -f (y^ - 77)2] 


/ 


, =^(.> = 2 I 0l 

0x 0y 


^"‘ (Xi - ^)2 + (y^ - r?)2j 


ds 


/ 


^ , (Xj - l)(y. - 7?) 

K. =-i^(d..) = -2 ■ 

' Ui - {)" -^ (Vi - »r 


ds 


The evaluation of these integrals is given in reference 9. 



APPENDIX B 


NUMERICAL REPRESENTATION OF THE FUNCTION OF PLASTIC STRAINS g(x, y) 


Before we can proceed with the numerical solution, it is necessary to represent the 
function of plastic strains g(x, y), given by equation (2) or (3), by corresponding finite- 
difference equations. 

The finite-difference net for grid station (r, s) is shown in figure 28. For a given 
function f = f(x, y), by use of central differences, we can obtain the following expres- 
sions for derivatives of this function; 


fx = e,-l, s('r-l, s - 'r, s ‘ ‘r, s' 


^y ^r, s-l^^r, s-1 ~ ^r, s^ ^r, s+l^^r, s+1 ~ ^r, s^ 


^xx ‘^r-l, s^^r-1, s ” ^r, s^ “r+1, s^^r+1, s " ^r, s^ 


Vy s-l^^r, s-1 ~ ^r, s^ ®r, s+l^^r, s+1 ” ^r, s^ 


MBl) 


^xy ^r-1, s-l^r-1, s-1 ^r-1, s^r-1, s ^r-1, s+l^r-1, s+1 


^r, s-l^r, s-1 ^r, s^r, s ^r, s+l^r, s+1 ^r+1, s-l^r+1, s-1 


^r+1, s^r+1, s ^r+1, s+l^r+1, s+1 


where subscripts x, y denote differentiation with respect to variable x or y, r is a 
row index, s is a column index, and the coefficients are given by the following relations; 

' I 
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s 


Ax Ax„ 


*^r+l, s 


Ax AXt 


%s-l 


Ay Ayj^ 


' Ay Ay, 


^r-1, s 


Ax AXr 




r, s-1 


Ay Ay^ 


r, s+1 


AyAyj^ 


^r-1, s-1 “ ^r-1, s^r, s-1 
^r-1, s+1 “ ^r-1, s^r, s+1 
^r+1, s-1 “ ^r, s-l^r+1, s 
^r+1, s+1 “ ^r, s+l^r+1, s 
^r-1, s ~ ”^r-l, s-1 ' ^r-1, s+1 
^r, s-1 “ ~^r-l, s-1 ' ^r+1, s-1 


r, s+1 "^r-1, s+1 ” ^r+1, s+1 


^r+1, s ~^r+l, s-1 " ^r+1, s+1 


^r, s ^r-1, s-1 ^ ^r-1, s+1 ^r+1, s-1 ^r+1, s+1 
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where Ax, Ay, Ax,p, AXg, Ay^^, and Ayj^ are distances as defined in figure 28. When 
relations (Bl) and (B2) are used, the function g(x, y) can be expressed by following 
finite -difference expressions: 


g(x, y) = L.,_ ^(€P t AcP) - ^ 

I - n'^ L 


+ Ae^ 
X ^x 


'^r+l, s(^x 






j ^ ^ "r, s-l(^x ^ g_l ■ s-1 ®r, s+l^(^x 


®r,s+l(^x'^^^x)^ + "r-l,s(^y 

- (“r-l, X + “ryl, s>('? + 3 + “ryl. sf ? * 3 ] 

* s-l('xy + 3_1_ 5./ ’'r- 1 , s(‘xy '^^y)r-l, 

1 ” /i- 


■"r- 1 , s +1 f 


^^xv ^^xv) ■*■ s-lf^xv ■*■ ^^xv) 

^y^r-l,s+l >^’Si\xy 


^r, s(^xy s ^ s+l(^xy ■*■ 

■*■ ^r+1, s-l(^xy + g_l ^r+1, s(^xy + g 

^ ^r+1, s+l(^xy J 


Plane strain 
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g(x, y) = -E 


("r, s-1 ('x * - (“r. s-1 * “r, s+l>(^x + ^ 


- <“r-l, s + “r+1, s>('? * . * “r+1, b( 4 * ^‘y) 


r+1, s 


■ ^ _^r-l,s-l(^xy ■*''^^xy)j._j^g_j ^r-1, s(^xy + 


-1, s+l(" 


+ yr_i ^.ifeS„ + AeP 




feP + AeP 


^r, s-1 


> (B4) 


^y/r-i,s+i 

s(^xy x,oTxy^ xxj, -J/r, s+1 

’'r+l, s-lf'^xy s-1 ^ ® ^^^y ^ ^^^y^r+l, 




P ) 

xyA 


+ ^r+l,s+l6xy*''«xy),^l_^^lj 


Plane stress 
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TABLE I. - ORDER OF STRESS SINGULARITY n AT THE TIP OF THE NOTCH 


FOR A SPECIMEN WITH A SINGLE EDGE NOTCH SUBJECTED TO 
PURE BENDING - PLANE STRAIN 
[Poisson's ratio ii = 0. 33.] 


Dimension- 
less notch 
depth, 
a 

Notch 

angle, 

a, 

deg 

Strain hard- 
ening pa- 
rameter, 
m 

Elastic 

Dimensionless load, q 

0. 4 

0. 5 

0.6 

0.7 

0.8 

0.9 

0. 3 

3 

0. 10 

0. 4999 


0. 488 

0. 490 

0. 487 

0.473 

0. 475 

. 3 

10 

. 10 




. 496 

. 497 

. 492 

. 480 

. 487 

. 5 

10 




0. 499 

. 496 

.480 

. 472 



. 5 

10 

. 10 



.496 

. 498 

. 480 

. 478 




TABLE n. - ORDER OF STRESS SINGULARITY n 
AT THE TIP OF THE NOTCH FOR A SPECIMEN 
WITH A SINGLE EDGE NOTCH SUBJECTED 
TO PURE BENDING - PLANE STRESS 


[Dimensionless notch depth a = 0. 3; notch angle 
a= 10°; strain hardening parameter m = 0. 10; 
Poisson's ratio p = 0. 33.] 


Elastic 

Dimensionless load, q 

0. 5 

0.6 

0.7 

0.8 

0.9 

0. 4999 

0. 419 

0. 434 

0. 448 

0. 451 

0. 458 





























Figure 5. - Distribution of boundary subdivisions and typical interior grid for elastoplasfic 
problem. (All distances are made dimensionless with respect to beam width w. ) 
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Figure 7. - Growth of plastic zone size with load for specimen with single edge notch subjected to pure 
bending. Plane strain; dimensionless notch depth a»0. 5; notch angle a = lOP; strain hardening 
parameter m -0.05; Poisson's ratio p-0. 33. 
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Figure 8. - Growth of plastic zone size with load in vicinity of notch for specimen 
with single edge notch subjected to pure bending. Plane strain; dimensionless 
notch d^th a = 0. 5; hotch angle a ■ 10°; strain hardening parameter 
m - 0.05; Poisson's ratio p " 0.33. 


1 



Dimenston 


Figure 9. - Growth of plastic zone size with load for specimen with single edge notch subjected to pure 
bending. Plane strain; dimensionless notch depth a - 0. 5; notch angle a - 10°; strain hardening 
parameter m -0.10; Poisson's ratio p = 0. 33. 




Dimension- 
load, 
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Figure 10. - Growth of plastic zone size with load in vicinity of notch for specimen 
with single edge notch subjected to pure bending. Plane strain; dimensionless 
notch depth a - 0. 5; notch angle a ■ 10°; strain hardening parameter 
m =■ 0. 10; Poisson's ratio p ■ 0. 33. 
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Figure 11. - Growth of plastic zone size with load for specimen with single edge notch subjected to pure 
bending. Plane strain; dimensionless notch depth a - 0. 3; notch angle o • 10°; strain hardening 
parameter m*0. 10; Poisson's ratio pi*0. 33. 
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Figure 12. - Growth of plastic zone size with load m vicinity of notch for specimen 
with single edge notch subjected to pure bending. Piane strain; dimensionless 
notch depth a "0.3; notch angle a - IcP; strain hardening parameter 
m » 0. 10; Poisson's ratio p - 0. 33. 






Dimension 


Figure 13. - Growth of plastic zone size with load for specimen with single edge notch subjected to pure 
bending. Plane stress; dimensionless notch depth ? ■ 0. 3; notch angle a ■ 10°; strain hardening 
parameter m - 0. 10; Poisson's ratio pi • 0.33. 
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Figure 14. - Growth of plastic zone size with load in vicinity of notch for a specimen 
with single edge notch subjected to pure bending. Plane stress; dimensionless 
notch depth a • 0. 3; notch angle a • lOP; strain hardening parameter m = 0. 10; 
Poisson's ratio u - 0. 33. 
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Figure 15. - Dimensionless equivalent stress contours in vicinity of notch for specimen 
with single edge notch subjected to pure bending. Plane strain,- dimensionless load 
tj - 0. 9; dimensionless notch depth a » 0. 3; notch angle a = 10°; strain hardening 
parameter m = 0. 10; Poisson's ratio p - 0. 33. 
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Figure 16. - Dimensionless equivalent stress contours in vicinity of notch for specimen 
with single edge notch subjected to pure bending. Plane stress; dimensionless load 
q=0.9; dimensionless notch depth a -0.3; notch angle a ■ 10°; strain hardening 
parameter m = 0. 10; Poisson's ratio p - 0.33. 







Figure 18. - Dimensionless x-directional stress distribu- 
tion along x-axis for specimen with single edge notch 
subjected to pure bending. Plane strain; dimensionless 
notch depth a - 0.5; notch angle a = 10°; strain hard- 
ening parameter m *0.05; Poisson's ratio p»0. 33. 



Figure 19. - Dimensionless y-directional stress distribu- 
tion along x-axis for specimen with single edge notch 
subjected to pure bending. Plane strain; dimensionless 
notch depth a -0.5; notch angle o = ICP; strain hard- 
ening parameter m -0.05; Poisson's ratio p - 0.33. 
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Figure 23. - Variation of dimensionless generalized stress intensity factor with load for 
specimen with single edge notch subjected to pure bending. Plane strain,- dimensionless 
notch depth a =■ 0. 5; notch angle a = 10°; Poisson's ratio p = 0. 33. 



Dimensionless load, q 

Figure 24. - Dimensionless piqne strain y-directional notch opening displacement for 
specimen with single edge notch subjected to pure bending. Dimensionless notch 
depth a = 0. 5; notch angle a = 10°; Poisson's ratio p = 0. 33; stress-strain curve 
given by figure 25. 
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Figure 25. - Stress-strain curve for aluminum 5083-0 used in test (private communication from 
R. T. Bubsey and M. H. Jones of NASA Lewis Research Center). Young's modulus of elas- 
ticity E ■ 7.79x10^ newtons per square centimeter (11. 3x# Ib/in. Poisson's ratio 


Dimensionless load, q 


Figure 26. - Dimensionless plane strain J integral for specimen with single edge notch 
subjected to pure bending. Dimensionless notch depth a ■ 0. 5; notch angle o = 10°; 
strain hardening parameter m '■ 0. 05; Poisson's ratio p - 0. 33. 
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Figure 27. - Variation of ratio of dimensionless J integral to square of dimensionless 
generalized stress intensity factor with load for specimen with single edge notch 
subjected to pure bending. Plane strain; dimensionless notch depth a -0.5; notch 
angle a - 10°; strain hardening parameter m-0.05; Poisson's ratio p = 0. 33. 



Figure 28. - Finite difference net for station (r, s). 
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